/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;

namespace dnAnalytics.LinearAlgebra.Solvers.Preconditioners
{
    /// <summary>
    /// An incomplete, level 0, LU factorization preconditioner.
    /// </summary>
    /// <remarks>
    /// The ILU(0) algorithm was taken from: <br/>
    /// Iterative methods for sparse linear systems <br/>
    /// Yousef Saad <br/>
    /// Algorithm is described in Chapter 10, section 10.3.2, page 275 <br/>
    /// </remarks>
    /// <include file='../../../../examples/LinearAlgebra/Solvers/Preconditioners/Diagonal.xml' path='example'/> 
    /// <include file='../../../../examples/LinearAlgebra/Solvers/Preconditioners/IncompleLu.xml' path='example'/> 
    public sealed class IncompleteLU : IPreConditioner
    {
        /// <summary>
        /// The matrix holding the lower (L) and upper (U) matrices. The
        /// decomposition matrices are combined to reduce storage.
        /// </summary>
        private SparseMatrix m_LuDecomposition;

        /// <summary>
        /// Returns the upper triagonal matrix that was created during the LU decomposition.
        /// </summary>
        /// <returns>A new matrix containing the upper triagonal elements.</returns>
        internal Matrix UpperTriangle()
        {
            SparseMatrix result = new SparseMatrix(m_LuDecomposition.Rows);
            for (int i = 0; i < m_LuDecomposition.Rows; i++)
            {
                for (int j = i; j < m_LuDecomposition.Columns; j++)
                {
                    result[i, j] = m_LuDecomposition[i, j];
                }
            }

            return result;
        }

        /// <summary>
        /// Returns the lower triagonal matrix that was created during the LU decomposition.
        /// </summary>
        /// <returns>A new matrix containing the lower triagonal elements.</returns>
        internal Matrix LowerTriangle()
        {
            SparseMatrix result = new SparseMatrix(m_LuDecomposition.Rows);
            for (int i = 0; i < m_LuDecomposition.Rows; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    if (i == j)
                    {
                        result[i, j] = 1.0;
                    }
                    else 
                    {
                        result[i, j] = m_LuDecomposition[i, j];
                    }
                }
            }

            return result;
        }

        /// <summary>
        /// Initializes the preconditioner and loads the internal data structures.
        /// </summary>
        /// <param name="matrix">
        /// The matrix upon which the preconditioner is based.
        /// </param>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="matrix"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="MatrixNotSquareException">
        /// If <paramref name="matrix"/> is not a square matrix.
        /// </exception>
        public void Initialize(Matrix matrix)
        {
            if (matrix == null)
            {
                throw new ArgumentNullException("matrix");
            }

            if (matrix.Rows != matrix.Columns)
            {
                throw new MatrixNotSquareException();
            }

            m_LuDecomposition = new SparseMatrix(matrix);

            // M == A
            // for i = 2, ... , n do
            //     for k = 1, .... , i - 1 do
            //         if (i,k) == NZ(Z) then
            //             compute z(i,k) = z(i,k) / z(k,k);
            //             for j = k + 1, ...., n do
            //                 if (i,j) == NZ(Z) then
            //                     compute z(i,j) = z(i,j) - z(i,k) * z(k,j)
            //                 end
            //             end
            //         end
            //     end
            // end
            for (int i = 0; i < m_LuDecomposition.Rows; i++)
            {
                for (int k = 0; k < i; k++)
                {
                    if (m_LuDecomposition[i, k] != 0.0)
                    {
                        double t = m_LuDecomposition[i, k] / m_LuDecomposition[k, k];
                        m_LuDecomposition[i, k] = t;
                        if (m_LuDecomposition[k, i] != 0.0)
                        {
                            m_LuDecomposition[i, i] = m_LuDecomposition[i, i] - t * m_LuDecomposition[k, i];
                        }
                        for (int j = k + 1; j < m_LuDecomposition.Rows; j++)
                        {
                            if (j == i)
                            {
                                continue;
                            }
                            if (m_LuDecomposition[i, j] != 0.0)
                            {
                                m_LuDecomposition[i, j] = m_LuDecomposition[i, j] - t * m_LuDecomposition[k, j];
                            }
                        }
                    }
                }
            }
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <returns>The left hand side vector.</returns>
        public Vector Approximate(Vector rhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (m_LuDecomposition == null)
            {
                throw new MissingMatrixException();
            }

            if (rhs.Count != m_LuDecomposition.Columns)
            {
                throw new NotConformableException();
            }

            Vector result = new DenseVector(rhs.Count);
            Approximate(rhs, result);
            return result;
        }

        /// <summary>
        /// Approximates the solution to the matrix equation <b>Ax = b</b>.
        /// </summary>
        /// <param name="rhs">The right hand side vector.</param>
        /// <param name="lhs">The left hand side vector. Also known as the result vector.</param>
        public void Approximate(Vector rhs, Vector lhs)
        {
            if (rhs == null)
            {
                throw new ArgumentNullException("rhs");
            }

            if (lhs == null)
            {
                throw new ArgumentNullException("lhs");
            }

            if (m_LuDecomposition == null)
            {
                throw new MissingMatrixException();
            }

            if ((lhs.Count != rhs.Count) || (lhs.Count != m_LuDecomposition.Rows))
            {
                throw new NotConformableException();
            }

            /*
            Solve:
            Lz = y
            Which gives
            for (int i = 1; i < matrix.RowLength; i++)
            {
                z_i = l_ii^-1 * (y_i - SUM_(j<i) l_ij * z_j)
            }
            NOTE: l_ii should be 1 because u_ii has to be the value
            */
            Vector rowValues = new DenseVector(m_LuDecomposition.Rows);
            for (int i = 0; i < m_LuDecomposition.Rows; i++)
            {
                // Clear the rowValues 
                rowValues.Clear();
                m_LuDecomposition.GetRow(i, rowValues);

                double sum = 0.0;
                for (int j = 0; j < i; j++)
                {
                    sum += rowValues[j] * lhs[j];
                }

                lhs[i] = (rhs[i] - sum);
            }

            /*    
            Solve:
            Ux = z
            Which gives
            for (int i = matrix.RowLength - 1; i > -1; i--)
            {
                x_i = u_ii^-1 * (z_i - SUM_(j > i) u_ij * x_j)
            }
            */
            for (int i = m_LuDecomposition.Rows - 1; i > -1; i--)
            {
                m_LuDecomposition.GetRow(i, rowValues);

                double sum = 0.0;
                for (int j = m_LuDecomposition.Rows - 1; j > i; j--)
                {
                    sum += rowValues[j] * lhs[j];
                }

                lhs[i] = 1 / rowValues[i] * (lhs[i] - sum);
            }
        }
    }
}
